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ABSTRACT 

Context. The exploitation of clusters of galaxies as cosmological probes relies on accurate measurements of their total gravitating 
mass. X-ray observations provide a powerful means of probing the total mass distribution in galaxy clusters, but might be affected 
by observational biases and rely on simplistic assumptions originating from our limited understanding of the intracluster medium 
physics. 

Aims. This paper is aimed at elucidating the reliability of X-ray total mass estimates in clusters of galaxies by properly disentangling 
various biases of both observational and physical origin. 

Methods. We use N-body/SPH simulation of a large sample of ~ 100 galaxy clusters and investigate total mass biases by comparing 
the mass reconstructed adopting an observational-like approach with the true mass in the simulations. X-ray surface brightness and 
temperature profiles extracted from the simulations are fitted with different models and adopting different radial fitting ranges in order 
to investigate modeling and extrapolation biases. Different theoretical definitions of gas temperature are used to investigate the effect 
of spectroscopic temperatures and a power ratio analysis of the surface brightness maps allows us to assess the dependence of the 
mass bias on cluster dynamical state. Moreover, we perform a study on the reliability of hydrostatic and hydrodynamical equilibrium 
mass estimates using the full three-dimensional information in the simulation. 

Results. A model with a low degree of sophistication such as the polytropic /3-model can introduce, in comparison with a more 
adequate model, an additional mass underestimate of the order of ~ 10% at rjoo and ~ 15% at r2m- Underestimates due to extrapolation 
alone are at most of the order of ~ 10% on average, but can be as large as ~ 50% for individual objects. Masses are on average biased 
lower for disturbed clusters than for relaxed ones and the scatter of the bias rapidly increases with increasingly disturbed dynamical 
state. The bias originating from spectroscopic temperatures alone is of the order of 10% at all radii for the whole numerical sample, 
but strongly depends on both dynamical state and cluster mass. From the full three dimensional information in the simulations we find 
that the hydrostatic equilibrium assumption yields masses underestimated by ~ 10 - 15% and that masses computed by means of the 
hydrodynamical estimator are unbiased. Finally, we show that there is excellent agreement between our findings, results from similar 
analyses based on both Eulerian and Lagrangian simulations, and recent observational work based on the comparison between X-ray 
and gravitational tensing mass estimates. 

Key words. Galaxies: clusters: general - X-rays: galaxies: clusters - Methods: numerical 



1. Introduction 

Galaxy clusters are the largest virialized structure known in 
the universe. According to the hierarchical clustering model of 
structure formation they form by the gravitational collapse of the 
rare peaks of the primeval density field, on scales of the order of 
~ 10 Mpc. Within this scenario their formation and evolution is 
a sensitive function of the cosmological matter density parame- 
ter fini and the mass fluctuation amplitude erg, where erg is the 
rms linear fluctuation on scales 8 /j"' Mpc. 

Measurements of their evolution rate can be used to asses the 
growth in mass of such structures, thereby providing a powerful 
method t o constrain the geometry and matter content of the uni- 
verse (see lVoitll2005L and references therein). Moreover, because 
of the spatial extent of the collapse scale, the cluster baryonic 
fraction /b is expected to be close to the cosmic value Qb/Qm- 
Measurements of /b at high redshifts can be used to derive con- 
straints on the equation of state of the dark energy dHaiman et alj 
l200lt iMaiumdar & MohilllOOl lAllen et aljllooj . 

The importance of cluster of galaxies as cosmological probes 
will be further strengthened with the upcoming high redshift sur- 
veys. This is of particular relevance in the new era of precision 
cosmology, in which studies of cluster evolution will provide 



independent tests with which to compare constraints on cosmo- 
logical models extracte d from observations of the cosmic back- 
ground radiation (e.g. ISpergel et al.l l2007l) and distance mea- 
surements of high reds hift supernovae (e.g. iTonrv et akl 120031 : 
iRiess et aLllIOOllIOOTl) . 

From the scenario here outlined it follows that in or- 
der to use cluster of galaxies as cosmological probes it is 
crucial to accurately measure their baryonic and total grav- 
itating mass. The methods used to derive cluster masses 
are mainly based on th e velocity dispersion of the optical 
galaxy populations (e.g. iBiviano & Girardil 120031 : iRines et akl 
120031) . observ ations of the X-ray emitt in g intracluster medium 
(ICM ) (e.g. iFinoguenovet al.f [20011: iReiprich & Bohringeii 



'2002'; 'Ettori et al.' '20021 lArnaud. Pointecouteau & Pratt! 12005 . 
Vikhlinin et al. 2006 ) and on gra vitational lensing (e.g. 
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iSmith et ail2(j02l:iMahdavi et alJl2008l) . 

Accurate mass estimates derived from X-ray data are based 
on the assumptions that both the total potential and the ICM dis- 
tribution are spherically symmetric and that the ICM is in hy- 
drostatic equilibrium in the cluster potential well. The latter as- 
sumption is usually justified by the fact that the estimated ICM 
soun d crossiiig tim es are short when compared against cluster 
ages (ISarazinlll986l) . 
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Under these assumptions the ICM is a faithful tracer of 
the underlying matter distribution and the total mass profile 
can be determined from the gas radial density and tempera- 
ture profiles. The density profile is recovered from deprojecting 
of the surface brightness profile, as measured from X-ray flux 
maps, whereas knowledge of the temperature profile requires the 
availability of spatially resolved spectroscopy. High quality data 
taken by the present generation of X-ray telescopes, Chandra 
and XMM-Newton, allows for nearby clusters accurate measure- 
ments of these quantities out to a large fraction of the clus- 
ter virial radiu s jMarkevit ch et al. 1998; Vikhlini n et al.l 
Piff'aretti et al. 2005; Sand erson et all 120061: IPratt et al 
iLeccardi & Molendi 2008) . 



2005 



2007 



The reUabiUty of cluster mass estimated through the X-ray 
method can be accurately studied using N-body/hydrodynamical 
simulations. The principal benefit over analytical methods is that 
the gas evolution can be treated self-consistently. The validity 
of the numerical approach is supported by X-ray observations, 
which show the existe nce of complex thermal struct ures and 
of merging activity (see lMarkevitch & Vikhlininll2007L for a re- 
view). 

Since early pioneering studies, hydrodynamical simulations 
have become a widely used tool to investigate cluster for- 
mation an d evolution in different cosmological scenarios (cf. 

I2005h . The numerical resolution of the simulations and 
the modelization of the cluster gas physics has been improved 
over the years. The latter now incorporates radiative cooling 
jYoshikawa. Jing & Sutoll2000|; iLewis et al.1 120001: iPearce et all 
120001: iMuanwong et al J 120021: loave et al.ll2002l) . metal_enricli 



ment o f the ICM by supernov ae and energy feedback ( K av et al 
20031: lTornatoreet"alT 120031: IValdarnini 2003; Borgani et al 



200llKav et al.ll2004l:lKravtsov et al.ll2005[ 



The accuracy of cluster X-ray mass estimators has been 
tested by mean s of N-body/hydrodynamical simulations in a va- 
riety of pap ers ('Evrard''l990'; Evrard et al. 1996;'Kav et al.''2004); 
iRasia et al.1 12006,; ,Kav et al.,,2007^ .Nagai et al., .2007a) . .Evrardi 
d 19901) first pointed out the existence of a significant bias in 
the binding mass es timates when using the isothermal j6-model. 
lEvrard et alj ( 1 19961) confirmed that the source of this discrep- 
ancy is related to the isothermal and hydrostatic assumptions. 
The lack of validity of the hydrostatic assumption is observa- 
tionally motivated by optical and X-ray maps, which show the 
existence of substructure with an ongoing merger ac tivity, and 
is numerically supported by a number of aii thors (e.g. lKav et al.l 
120041; iRasia et alll2004l; iNagai et aljr2007ah . who found that in 
the simulations the ICM is not perfectly in hydrostatic equilib- 
rium. This implies the presence of residual gas bulk (laminar) 
and turbulent (random kinetic) motions, and leads to an underes- 
timate of the masses because of additional non-thermal pressure 
support which is not accounted for by the hydrostatic equilib- 
rium equation. 

With respect the isothermal j6-model the modeliza- 
tion of the ICM h as been significantly improved (e.g. 
IVikhlinin et alj l2006h with observational progresses, which 
showed the existence of temperature profiles declining with 
radius (iDe G randi & Moiend? '200?, 'Vikhlinin et al.' ^200l 
iPiffaretti et ai . 2005; Sanderson et al. 2006; Pratt et al. 2007). 
These features are well reproduced out of the core radii in hydro- 
dynamical simulations which incorporate cooling and feedback 



dynamical simulations which mcorporate cooling and leedbacK 
jMuanwong et al 1 120021; iKav et al.ll200l iTornatore et al.1 120031: 
IValdarninill2003l;lBorgani et alj|2004~ 



In order to properly assess the reliability of cluster X-ray 
mass estimators it is however necessary to construct mock ob- 
servations of simulated clusters which must reproduce spectro- 



scopic measurements as expected from X-ray telescopes. This 
is motivated by the presence of complex thermal structures 
in the ICM, which bias the (measured) spectral fit tempera- 
tures towards lower values than the averag e emission weighted 
cluster temperatures defined theoretically ( Mathiesen & Evrard 
I2OOII; iMazzotta et al.l [20041; IValdarninill2006l) . The dependence 
of X-ray mass estimators on spectral biases and other sys- 
tematics has been investigate d through hydrodynami cal simula- 
tions by a number of authors (Rasia et al. 2006 : Kav et al. 2QQ% 
INagai et al.ir2007al; IJeltema et al.ll2008h . lNagai et all(i2007a) ar- 
gued that mass estimates are biased low (~ 5 - 20%) even for 
clusters identified as relaxed. 

In order to properly disentangle observational biases which 
arise from spectroscopic measurements from those due to in- 
complete relaxation of the gas or from the ones caused by an in- 
accurate modeling of the radial profiles, it is however necessary 
to derive X-ray mass estimates from a large sample of simulated 
clusters. This is the main goal of this paper, in which we apply 
different X-ray mass estimators to a large set of high-resolution 
hydrodynamical simulations of galaxy clusters. The physics of 
the gas includes radiative cooling, star formation, chemical en- 
richment and energy feedback. The sample comprises ~ 100 
clusters, the size of the sample being a critical quantity in or- 
der to extract sub-samples large enough to give a meaningful 
statistics. 

We discuss the dependence of the mass bias at different radii 
upon the adopted analytical models and the chosen radial range 
used to perform the fits of the profiles, the cluster dynamical state 
as well as the impact on the mass bias which follows from the 
use of spectral temperatures. As a statistical indicator to quantify 
the cluster dynamical state through the analysis of X-ray m aps 
we adopt the power ratio method (see IJeltema et al.ll2005l and 
references therein). We also investigate the limit of applicability 
of the dynamical equilibrium equation when used to recover the 
cluster true masses in the presence of significant non-thermal gas 
pressure. 

The paper is organized as follows. In Sect. |2| we present the 
procedure for simulating the cluster sample. In Sect. [3] we de- 
scribe how we generate and analyze the synthetic X-ray obser- 
vations that are used in Sect. |4| to recover the total mass dis- 
tribution. In Sect. |5] we investigate the reliability of the hydro- 
static and hydrodynamical mass estimators from the full three- 
dimensional information provided by the simulations. Finally, 
we discuss our main results and present our conclusions in Sect. 



2. N-body/SPH simulations 

The considered cosmological model assumes a flat CDM uni- 
verse, with matter density parameter Q.^ - 0.3, Qa = 0.7, 
Qb - 0.0486 and h = 0.7 is the value of the Hubble constant in 
units of 100 km sec^'Mpc"'. The power spectrum has been nor- 
malized to cTg - 0.9 on a 8/1"' Mpc scale at the present epoch 
and the primeval spectral index n is set to 1 . 

The simulation ensemble of galaxy clusters is constructed 
according to the following procedures. A low-resolution N-body 
run is first performed starting from an initial reshift Zi = 10 
in a box of comoving size L, using a P3M code with A^p par- 
ticles. Clusters of galaxies are identified at z = using a friends- 
of-friends (FoF) algorithm. The virial mass and radius are re- 
lated by Mvh- = (47r/3)QmPcAcr^;^., where Ac - 187 O^^-^^ 
for a flat cosmology and pc is the critical density. This is de- 
fined as pe(z) = 3H{zf/87TG, where H{zf = H?.E(zf and 
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E{z)^ - fimCl + zy" + Qa- In general, the fiducial radius is de- 
fined such that the average of the total density within that radius 
is A times the critical density, i.e. M(< ta) = 47rr^ Apc(z)/3. 

The simulation ensemble is constructed by combining two 
distinct samples SI and S2. Table[T]lists the most relevant simu- 
lation parameters of the two samples SI and S2. 

Table 1. Main parameters of the two cluster samples SI and S2: 
comoving box size L of the cosmological A^-body runs, number 
of corresponding particles Np, number Nd of the most massive 
clusters extracted from the sample at z = 0, and virial mass M,,,,- 
of the most (least) massive cluster of sample S1(S2). 



Sample 


L 






M,„, 




[/r'Mpc] 






[io"'r'Mo] 


SI 


400 


168^ 


33 


15 


S2 


200 


84^ 


120 


1.5 



A sample is constructed by extracting at z = the Nd most 
massive clusters which are found in a cosmological run of size L 
and A^^ particles. For sample SI the box size and the number of 
particles are twice and eight times those of sample S2, respec- 
tively. The value of for sample S 1 is chosen such that the last 
cluster of the sample has its virial mass above that of the first 
cluster of sample S2. The mass range for the virial masses of the 
combined samples spans a decade. The random realization of the 
initial density perturbations are different in the two cosmological 
simulations. 

After competition of the cluster selection, the clusters of the 
ensemble are resimulated individually using high-resolution hy- 
drodynamic simulations in physical coordinates. The initial con- 
ditions of each hydrodynamic simulation are set as follows. The 
particles of a cluster lying at z = within a distance r^^ from the 
cluster center are tagged and located back at an initial redshift 

1 /3 

Zin = 49 and a cubic region of size ^ 25 - 50 Mpc oc MJ^^ , 
which contains these particles, is placed at the cluster center. The 
original Fourier modes of the cosmological simulation are then 
used to perturb the positions of a uniform lattice of Nl - 51^ gas 
particles set inside the cube and high frequency waves are added 
to the original modes to sample the new Nyquist frequency. The 
positions of dark matter particles are found similarly. Particles 
with perturbed positions which lie within a inner sphere of ra- 
dius Lc /2 from the cluster centre are kept for the hydrodynamical 
simulations. For these particles the masses are set in proportion 
to Qb and (fim - Qb) for gas and dark matter particles, respec- 
tively. Tidal fields are modeled by adding to the inner particles an 
external shell of low-resolution dark matter particles. The shell 
has an outer radius and the mass of a particle is 8 times the 
sum of the masses of a gas and dark matter particle of the inner 
region. 

The hydrodynamic simulations are run using a multistep 
TREES PH code in which the gas entropy is explicitly con- 
served jGoodman & Hernquistlll99l]) . The simulations contain 
^ 70, 000 gas and dark matter particles in the inner region and a 
similar value of low -resolution dark matter particles in the outer 
shell. The mass of the gas particles ranges from mg ^ 5 x IO^Mq 
for the most massive cluster of the ensemble, down to ^ 
6 X IO^Mq for the least massive cluster This mass resolution can 
be considered adequate for the present purposes, as suggested 
from analyzing th e stability of gas profiles of simulated clusters 
d Valdarninill2002[) . The gravitational softening parameter of gas 
particles is set to Sg - 25 kpc and 15 kpc, for clusters of sample 



S 1 and S2, respectively. For the dark matter particles the soften- 

1 /3 

ing is rescaled according to Si oc mj ' , where mi is the mass of 
the particle /. The softenings are comoving out to z = 20, after 
which are kept fixed in physical coordinates. 

The physical modeling of the gas includes radiative cooling, 
which depends on the gas temperature and metallicity. Cold gas 
in high density regions is subject to star formation and gas parti- 
cles are eligible to form star particles. At each timestep gas parti- 
cles neighboring a star particle are heated by supernova (SN) ex- 
plosions of type II and la. The gas is also metal enriched through 
SN explosions. The energy and metal feedback are calculated 
according to the stellar lifetime and initial mass func tion. A de- 
tailed description of the feedback recipes is given in I Valdarninil 
(l2003h . The hydrodynamic variables of a simulated cluster are 
stored at run time at various redshifts in the interval from z = 2 
down to z = 0. 



3. Simulation and analysis of X-ray observations 

At z = we extract a total of 153 clusters. For this sample we 
compute spectroscopic-like global temperatures and select ob- 
jects with T^\{< rsoo) > 2 keV (see Sect. l3.4l below for the defini- 
tion of Tsi). This temperature selection is adopted because when 
generating spectrosco pic temperature profile s we use the ap- 
proximation derived bv lMazzotta et al.l (|2004|) . which was devel- 
oped for continuum-dominated spectra and is therefore not ac- 
curate for low temperature systems. Our final sample comprises 
~ 100 temperature-selected clusters spanning a mass range of 
8.2 X W^H'^Mq < M{< r2oo) < 1.2 x IO^^/j-'Mo. 

For each simulated cluster we generate three independent ob- 
servations by considering three orthogonal projections, thus ob- 
taining a total of ~ 300 observed objects. In following sections 
we describe how we generate and analyze these mock X-ray ob- 
servations. Our an alysis may be vie wed as complem e ntary t o 
those presented in iRasia et al] ( |2006[) and iNagai et al.l ( l2007al) . 
The main difference being the large size of the sample used here, 
which allows to draw results of much higher statistical signifi- 
cance. 

3.1. X-ray maps 

Simulated surface brightness maps are obtained by first choosing 
a line of sight and then locating the origin at the cluster center in 
the plane orthogonal to the line of sight, the latter is defined as 
the location where the gas density reaches its maximum value. 
Throughout the paper we assume that this position also coincides 
with the peak of the X-ray emission. 

For a given line of sight and redshift z the X-ray surface 
brightness is defined as: 

Sx(x,y)^ —^—j r sxdl, (1) 
47r(l +z)4 J_co 

where the integral is along the line of sight, x and y are cartesian 
coordinates on the chosen plane and ex is the X-ray emissivity. 
The latter quantity can be calculated as ex = Apg(x), which dif- 
fers from the true emissivity aside from a constant factor. Here 
Pg(x) is the gas density, A = A(T,Z,x, Ei, E2) is the cooling 
function, and T and Z are the gas temperature and metallicity, 
respectively. The quantities £1 and E2 define the energy band 
[E1-E2] used in the X-ray flux measurement. Because of the 
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Lagrangian nature of SPH simulations, the emissivity ex(x) is 
expressed as a summation over particles: 

ex(x) = y —K{Tj, Zj, X, El , E2)pI /x)H'(x - xj, hj) , (2) 

j Pg.j 

here the subscript j denotes the values of the quantity at the posi- 
tion of the particle j, W is the SPH smoothing kernel, and hj the 
smoothing length of particle j. From this expression for ex(x) it 
is possible to evaluate the X-ray surface brightness Sx{x,y) by 
performing the integral ([1]) along the line of sight axis, z- For our 
purposes the integral ([U is replaced by 



1 



(3) 



where e^^ is defined according to eq. in which the kernel 
is the 2D Gaussian kernel W2D - exp[-{ yjx^ + / hjf] / nh'j ^ 
and/i^ ^ 0.82/!/. 

For each simulated object we generate three surface bright- 
ness maps by choosing three orthogonal projections. Motivated 
by observations we evaluate the maps in the soft energy band 
[0.5-2] keV. Since particles with temperatures below 10^ K do 
not contribute to the gas emissivity in this energy band, they are 
not taken into account in the computation. 

Cold and dense ICM clumps can affect total mass recon- 
structions because they produce pronounced irregular features 
in azimuthally averaged brightness and temperature profiles, and 
introduce a systematic bias in the spectral temperature determi- 
nation. We therefore remove small scale clumps from the con- 
structed flux images by performing a masking procedure as done 
with real X-ray observations. This is accomplished by generat- 
ing [0.5-2] keV maps on a 1024 x 1024 grid of size 2.1 x ^200- 
The pixel size of the maps is ~ 5 and 2 kpc for the most and 
least massive clusters in our sample, respectively. Pixel sizes are 
similar to those adopted in observations and sufficiently small 
as to allow a complete identification of all small-scale clumps 
resolved in our simulations. Clumps are d etected and remove d 
following the procedure implemented by Vik hlinin et al.l (Il998h : 
images are decomposed into wavelets of pixel scale 2", with 
n = 0, 1 , ... 6, with significance threshold set to 5 in units of the 
rms level. Particles which lie in those pixels which are part of a 
statistically significant structure on a given scale are tagged and 
identified as part of a clump. These particles are then removed 
from all the SPH summations. In three dimensions clumps are 
identified by those particles which are part of a two dimensional 
clump in any of the planes. The wavelets algorithm commonly 
detects also the central peak of the surface brightness distribu- 
tion. In order to avoid the unwanted masking of the cluster emis- 
sion in the central region we exclude the inner 0.2 X ?-2oo cir- 
cular aperture from our masking procedure. In addition, we do 
not impose any limiting flux fx for the detection of small scale 
clumps. We motivate this choice with the finding s that se tting 
fx to ~ 3 X 10^'^erg /cm'^ as in iNagai et"all ( 12007 ah . has 
fully negligible effects. Finally, we compute the amount of gas 
mass removed by the masking procedure and find that even in 
the most clumpy clusters it is less than a few percent of the total 
gas mass. If not stated explicitly, all the results presented in the 
following are derived by including the masking procedure. 

3.2. Power ratios 

Total mass determinations from X-ray observations can be heav- 
ily affected by the cluster dynamical state since the method 



relies on the assumption of hydrostatic equilibrium. The dy- 
namical state of clusters is related to the amount of substruc- 
ture present in their X-ray surface brigh tness distribution (e.g. 
iRichstone et"ai]|I992l:lEvrard et al.lll993h and various statistical 
mea sures have b een proposed to quantify cluster substructure 
(e.g. lBuot3l2002l and references therein). 

In this work we adopt the power ratio method dBuote & Tsail 

[T995I) as a statistical indicator of the cluster dynamical state, 
since it is widely us ed to study cluster X-ray morphologies (e.g. 
iBuote & TsailI996l) . This method is expect to be statistically sig- 
nificant when applied to a large cluster sample such as the one 
studied here. It is in fact unaffected only by mergers along the 
line of sight which rarely occur 

According to the power ratio method, the X-ray surface 
brightness map Sxip, f), where (p, (p) are the conventional polar 
coordinates, is the source term of the pseudo potential ^'{p, (/?) 
which satisfies the 2-D Poisson equation. The pseudo potential 
is expanded into plane harmonics and the m-th coefficients of the 
expansion are given by : 



d^x S x(x )R ' cos(mip ), 



/?,„ = ( d^x Sx{x )R " sm(tiup ), 

'<«<„, 



(4) 
(5) 



where x - (p,ip) and the integration is over a circular region 
with aperture radius Rap. The m-th power ratio is then defined as 



n„,{Rap) = logio(P,„/Po) , 



where 



1 



PnARap) = ^ (<+/?;;,) m > 0, 

Pq = [aoln(7;„/,/kpc)]2. 



(6) 

(7) 
(8) 



The power ratios UmiRap) are then indicators of the amount 
of structure present on the scale of the aperture radius /?ap- The 
values of depend on the choice of the coordinate system. 
For a fully relaxed configuration — > -00. As indicates 
asymmetric distributions, we adopt it here as a measure of the 
amount of substructure present in a cluster 

For each of the ~ 300 [0.5-2] keV surface brightness maps 
(see Sect. 13.11 ) we set the origin of coordinates at the peak of 
the X-ray emission. In addition to being fully consistent with 
the procedure employed in the computation of azimuthally av- 
eraged brightness and temperature profiles, this choice of coor- 
dinate system enables us to detect P3 values different from zero 
even in the case of bimodal clusters with nearly equal size com- 
ponents. For the same configuration all the odd moments would 
vanish if the coordinate system was the frequently adopted cen- 
troid fra me, which is the co ordinate system defined such that 
Pi = dBuote & Tsailll99"5l) . We utilize unmasked maps, since 
the goal of the power ratio analysis is to measure substructure, 
and for each one we compute n3(/?ap) in correspondence of the 
same radii at which we evaluate the mass biases (r25oo, '"500, and 
r2oo, see Sect. |4]i. The various integrals involved in the compu- 
tations are performed according to the SPH prescription. For the 
simulated clusters of our sample, the values of 113 we found lie in 
the range between ~ -(6 - 4) for a strongly asymmetric distribu- 
tion and ~ -(12 - 10) for a cluster with a relaxed configuration. 

In order to perform a statistical analysis of mass determi- 
nation biases as a function of substructure we extract four sub- 
samples in the following way. For a given overdensity, we first 
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construct the cumulative distribution of the sample values of 113. 
Then we identify the synthetically observed clusters in four sub- 
samples (113 classes, hereafter): first quartile (below 25%), sec- 
ond quartile (between 25% and 50%), third quartile (between 
50% and 75%), and forth quartile (above 75%). In the follow- 
ing these four sub-samples will be referred to as qi, q2, ^3, and 
q4 class, respectively. By construction the most relaxed clusters 
belong to the qi class and the most disturbed ones to the ^4 class. 

3.3. Computation an modeling of surface briglitness profiles 

The azimuthally averaged surface brightness profiles 5'x('") are 
computed from the [0.5-2] keV band surface brightness maps as 
follows. Each map is binned using a grid with circular geometry 
and coordinates (p, 0), the origin is set at the position of the X- 
ray emission peak. The grid points are uniformly spaced, linearly 
in the angular coordinate, and logarithmically in the radial coor- 
dinate. The range of the spatial coordinates is between 2 x 10""* 
and 1.5 in units of r20(). There are Np - 140, A^^ = 20 points in 
the coordinate intervals. The number and spacing of the radial 
grid points as chosen such that the surface brightness binning in 
the radial range used in the fits (see Sect. |4] below) is similar to 
the one adopted in the analysis of real X-ray data. The value of 
Sx(r) at a given projected radius r is given by averaging over 
the azimuthal values of the annulus. These values are also used 
to define a surface brightness dispersion, which is later used to 
define a weight in correspondence of 5'x('") when the profile is 
fitted. 

In order to explore the influence of the models adopted in 
the fitting of surface brightness radial profiles on total mass esti- 
mates, we adopt two different models. In the simplest and most 
widely adopted procedure used in the determination of the to- 
tal gravitational mass of clusters t he surface brightness profile 
is mo deled by a single j6-model dCavaliere & Fusco-Femianol 
fT976h . i.e. 



-3/3+0.5 



(9) 



The single y6-model is convenient because it is simple to depro- 
ject it and obtain the ICM density profile needed to estimate the 
total mass from the hydrostatic equilibrium. 

As a second parametrization we ado pt an extended ^-mode l 
in which the profiles are modeled as in I Vikhlinin et al ] (I2006h . 
The radial dependence of Pg(r) is given by the functional form: 



1 



«'°(1 + r2/r^yP-»n (1 + ry/r]r/y' 



(10) 



where the additional parameters with respect the standard 
jS-model allow a much more accurate modeling of the density 
slope variations. The component present in in Eq. (fTOl i which 
takes into account for a central surface brightness excess is ne- 
glected (i.e. we set a - 0), since, as discussed below, we are 
excising the very central regions in the fitting procedure. The 
model parameters are recovered by fitting the projection of Eq. 
([Tol l along the line of sight against the surface brightness pro- 
fil es in the [0.5-21 keV e nergy band. Following the suggestions 
of lVikhlinin et all (l2006l) we set in Eq. (fTOli y = 3. 

Notice that for both parametrizations the total mass derived 
from the hydrostatic equilibrium equation (Eq. (fl4l i below) does 
not depend on the central gas density value. 

As commonly found in simulated clusters, the gas properties 
in the central regions of our simulated objects are in disagree- 
ment with observations. The pronounced steepening of the gas 



density towards the center in low temperature systems and the 
sudden temperature drop in the very inner regions are in fact not 
observed in real clusters. We therefore exclude the inner region 
within riow = 0.1 x ^200 from the surface brightness fits. We em- 
phasize that the precise value for riow is not relevant, since our 
goal is to quantify the bias in the mass reconstruction at much 
larger radii. 

3.4. Computation an modeling of temperature profiles 

From the three-dimensional gas temperature distribution T^'^ix) 
measured in the simulated clusters we compute both three- and 
two-dimensional (projected) radial profiles. 

A three-dimensional radial temperature profile can be de- 
fined as 



(11) 



where is a weight function and the volume integral is 
over a spherical shell of thickness Ar located at distance r 
from the cluster center Common choices for W are the gas 
density (mass-weighted temperatures, W = Pg) and the X-ray 
emissivity (emission-weighted temperatures, W - sx). In this 
work we consider mass-weighted temperatures and the 
weight function W = p^T'^^"^, which defines the c orresp onding 
spectroscopic-like temperature r^i . iMazzotta et al] ( |2004|) found 
that in the continuum regime (T > 2keV) this choice of the 
weight function provides accurate approximations of spectro- 
scopic temperatures obtained from X-ray observations (i.e. de- 
rived by fitting an integrated spectrum with a single-temperature 
model). A more sophisticated method would be based on the 
temperature determinations from either projected spectra (in 
the case of projected temperature profiles) or from deprojected 
ones (in the case of deprojected temperature profiles) com- 
puted from the simulation outputs. Given the size of our sample 
and being this procedure computationally expensive, we adopt 
spectroscopic-like temperature profiles. 

For the cluster under consideration we evaluate the weighted 
temperature profiles by replacing the integrals in Eq. ( fTTT i. ac- 
cording to the SPH scheme, by a summation over gas parti- 
cles. The value of the radial temperature profile T^{r) in cor- 
respondence of each shell is defined by averaging over a set of 
(0, (j)) - 40x40 grid points uniformly spaced in angular coor- 
dinates. The radial coordinate of the shells is binned as in the 
computation of the surface brightness. For each object we there- 
fore obtain two types of three-dimensional temperature profiles: 
a mass-weighted profile and a spectroscopic-like one. 

The projected temperature profiles T^{r) are computed in 
the same way as the three-dimensional profiles, but by con- 
sidering a cylindrical geometry instead of a spherical one. 
Furthermore, only spectroscopic-like profiles are computed, 
since mass-weighted projected profiles are not meaningful in this 
case. For each cluster we compute three projected spectroscopic- 
like temperature profiles by consistently considering the same 
projections used in the computation of surface brightness pro- 
files. 

We compare projected to three-dimensional spectroscopic- 
like profiles and find that in the radial range [0.1-1] ^200 the dif- 
ference is generally small. The agreement between the two types 
of profiles steadily improves from disturbed to relaxed clusters 
and the difference is fully negligible for very relaxed clusters. 

In the computation of temperatures, in accord with the 
considered energy bands, we neglect in the SPH summations 
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particles with temperatures below 10^ K. For consistency, we 
adopt the same temperature cut-ofF in the computation of mass- 
weighted temperatures. Moreover, for the spectral weighting we 
restrict the summation to those gas particles for which kT > 0.5 
keV. 

It is crucial to remark that spectroscopic-like temperatures 
are good approximations of the spectroscopic temperatures as 
would be estimated from X-ray observations, whereas mass- 
weighted temperatures are a much more accurate indicator of 
the total binding mass, since they follow the virial relationship. 

As for the brightness profiles, we use two different models 
to fit temperature profiles and exclude the inner aperture within 
How =0.1 xr2oo from the fits. In the simplest case the gas temper- 
ature profile is modeled using the commonly adopted polytropic 
relation 



Tg oc p'g 



(12) 



with 1 < y < 5/3. 

In the second model we allow a larger parametric free dom 
and follow the modeling proposed bv lVikhlinin et alj(l2006l) . We 
model the temperature profile using 



Tg(_r) = Tgfl 



(13) 



which describes a power-law declining profile with a transition 
at r ~ ri. In the above model is absent the fcooi('") term, which ac- 
counts for the tempe rature decline in the ce ntral region because 
of radiative cooling ( IVikhlinin et alj ( l2006h ). This is motivated 
by our choice of the radial range over which the profiles are fit- 
ted. 

When modeling temperature profiles, a proper treatment 
would require to fit the model projection along the line of sight 
against the simulated spectral temperature profiles. Moreover, 
one has to take into account the presence of different tempera- 
ture components which can bias the single-temperature fit. Given 
the large size of our sample, we avoid the involved projection 
and recovery steps and we fit the model in Eq. ( fT3] l directly to 
the 3D spectroscopic-like temperature profiles. These are in fact 
expected to be in good agreement with deprojected temperature 
profiles derived from spatially resolved spectroscopy. However, 
in order to explore the crucial bias introduced by spectroscopic 
temperatures, we additionally apply the same procedure to 3D 
mass-weighted temperature profiles, since they provide the az- 
imuthal average of the true gas temperature distribution. 

4. Total mass determination from X-ray 
observations 

The total gravitating mass of a galaxy cluster is estimated from 
X-ray observations assuming that the ICM is in hydrostatic equi- 
librium and that its temperature and density distributions, as 
well as the total gravitational potential, are spherical symmetric. 
These assumptions lead to: 



M'=''(< r) 



kTgir) r 
Gjinip 



d In pg{r) d In Tg{r) 



d\nr 



d\nr 



(14) 



where M''^'(< r) is the estimated total gravitating mass within 
a cluster-centric distance r, Tg{r) and Pg{r) are the three- 
dimensional gas temperature and density profiles at the radius 
r, respectively, and G and nip are the gravitational constant and 
proton mass. We adopt the value ji - 0.58 for the mean molecu- 
lar weight of the gas. 



4. 1 . The mass bias 

In order to quantify the difference between estimated mass 
M'^^'(< r) and the actual mass of the simulated object M{< r), 
we define the mass bias as; 



M'^'iK r) - M(< r) 
M(< r) ■ 



(15) 



The true mass M{< r) is computed by summing the masses of 
all the particles (dark matter, gas, and star particles) within the 
radius r. 

In our observational-like analysis the functions Pg(r) and 
7'g(r) in Eq. ( fT4l i are modeled as described in Sects. [373] and l3. 41 
The various options for the choice of fitting functions and tem- 
perature profiles allow us to investigate and disentangle biases 
of different origin. In particular we explore three different cases; 

- Polytropic /3-model: for a given cluster and projection we 
fit the surface brightness profile and the corresponding pro- 
jected, spectroscopic-like temperature profile using Eqs. (|9]) 
and ( fT2l i. respectively. The derived mass bias is referred to as 

- Extended f3-model with spectroscopic temperatures: for a 
given cluster and projection we fit the surface brightness pro- 
file with Eq. ( fTOl i. The three-dimensional, spectroscopic-like 
temperature profile cluster is fitted using Eq. ( fT3] ). In this 
case the mass bias is denoted by feext/3T'°('")- 

- Extended fi-model with mass-weighted temperatures: same 
as the previous case, but the three-dimensional, mass- 
weighted temperature profile is adopted instead of the 
spectroscopic-like profile. The derived mass bias is denoted 

by ^exV3,T;;!i('-)- 

For each of the ~ 300 observations the derived mass bias is 
evaluated at three different overdensities, in correspondence of 
A — 2500, 500 and 200. Cluster-centric distances of r25oo and 
rsoo are usually probed by observations, with rasoo being well 
inside the largest accessible radius and rsoo being typically the 
largest distance where the ICM temperat ure can be reliably es 
timat e d through X-ray observations (e.g. De Gra ndi & Molendi 



Piffaretti 



A-ray 
etal.1 [ 



20051: iVikhUnin et all I2005e iPratt et al 



Snowden et all 120081: iLeccardi & Mole'ndTl2008l) . We 



12002: 
2007 

additionally choose rooa in order to extend our analysis at 
radii which will be accessible in the near future (e.g. see 
iReiprich et al.ll2068l) . to quantify the effects of extrapolation to 
large distances and, most important, to investigate the mass bias 
in dynamically different regions. 

We adopt three different values for the outer boundary of 
the region considered in the mass reconstruction by using data 
(i.e. surface brightness and temperature profiles) within r^^ - 
0.5, 0.75 and 1.05 x r2oo (the inner boundary is riow = 0.1 x r2oo, 
see Sects. [331 and |3.4t . The choice of rup is of paramount impor- 
tance to assess the dependence of the mass bias upon changes 
in the gas density slope, which can have a significant impact on 
mass measurement biases. The three different choices allow us 
to explore a limitation that in practice is imposed by the field of 
view of the detector and the background level of a given expo- 
sure. 

To summarize, for a given radial range ([0.1-0.5] r2oo, [0.1- 
0.75] r2oo, or [0.1-1.05] r2oo, with the latter denoted for sim- 
plicity by [0.1-1] r2oo hereafter), we fit surface brightness and 
temperature profiles with the analytic functions specific to the 
adopted model and evaluate Eqs. (O and ( [TSl ) at r25oo, '"500 > 
and r2oo- It is important to notice that by averaging over the 
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whole sample we find: < r25oo/'"200 >sampie= 0.31 and < 
'"50o/'"200 >sampie= 0.67. As a consequence M'^^'(< r), and there- 
fore b(r), is extrapolated using the best fit functions at rsoo and 
r2oo for the fitting range [0.1-0.5] r2oo and at r2oo for the fitting 
range [0.1-0.75] r2oo- No extrapolation is of course necessary 
when the fitting range [0.1-1] r2oo is adopted. 

The large size of the simulated mock sample (~ 300 obser- 
vations) allows us to use a statistical approach when exploring 
the dependence of the mass bias on various quantities. For any 
given mass determination method, overdensity at which the en- 
closed mass is estimated, radial fitting range, and YI3 class qi (or 
whole sample) the mass bias values are modeled using distribu- 
tion fitting. The data set is always modeled with both normal and 
WeibuU distributions in order to describe symmetric and asym- 
metric data distributions, respectively. We never find, however, 
strongly skewed distributions. Even when a Weibull distribution 
yields a better fit to the data than a normal distribution, its mean 
and standard deviation are only slightly different from those de- 
rived from the normal distribution best-fit. Since these small dif- 
ferences are completely irrelevant for our discussion, we quote 
only results from normal distribution fitting throughout the pa- 
per From the data set of values biir^,) (/ = /3, T~^; ext/3, T^^; 
ext/3, T^), distribution fitting yields the mean value m{bi(rA)) 
and standard deviation cr(bi(r/s,)). 

Notice that sample SI is extracted from a volume 8 times 
larger than the one for S2 (see Sect.|2]i. This implies that the clus- 
ter abundance measured in our sample is not correct. In order to 
check whether this issue has an important impact on our results 
we proceed as follows. We create a new sample made of sample 
S2 and 1/8 of randomly selected objects from sample SI, for 
which mean values and standard deviations of the mass bias are 
computed for the all the cases under consideration. The proce- 
dure is repeated many times and the derived values are compared 
with those computed for the original sample. In all the cases we 
find that the differences are fully negligible and we therefore re- 
port results for the original sample in the reminder of the paper 

Finally, it is important to notice that the four classes qi 
(see Sect. l3.2l ) are defined such that each class contains the same 
number of objects (~ 70). This guarantees a meaningful compar- 
ison between mean values and standard deviations derived from 
different 03 classes. 



4.2. Results 

Total cluster mass determinations from an X-ray observations, 
might be in practice affected by more that one source of bias. 
The number of different cases for which we derive total masses 
(fitting ranges, best fit functions, temperature definitions, etc.) 
are of course designed in order to investigate the various sources 
separately, but also to show the effect of their combination. Our 
results are reported in Table |2] and conveniently shown in three 
figures (Figs.[T][2] and[3j. Notice that in order to improve clarity 
the mean values (data points in the figures) and standard devi- 
ations (errorbars in the figures) of the various mass bias distri- 
butions are given in percent. The reported values show, in most 
of the cases, the combined effect of different biases and are very 
useful to estimate the total bias for a given mass reconstruction 
method and conditions. Nevertheless, some specific cases and 
the comparison of mass biases derived under different conditions 
allow us to quantify the effect of single biases separately. 
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Fig. 1. Summary of the results for the mass reconstruction using 
the polytropic j6-model and projected, spectroscopic-like tem- 
perature profiles. The three panels show quantities derived from 
modeling of the profiles in three different radial ranges. In each 
panel we show the mean (points) and standard deviation (error- 
bars) of the mass bias distribution at r2oo (stars), rsoo (diamonds), 
and r25oo (triangles). Open/filled symbols indicate quantities de- 
rived with/without extrapolation of the mass profile. The results 
from distribution fitting are shown for the whole sample (all) and 
when the four 113 classes are used as sub-samples. The results for 
different overdensities are shifted horizontally to improve clarity. 
Notice that the mass bias is shown here in percent. 

4.2.1. The modeling bias 

In order to disentangle the bias introduced by inaccurate mod- 
eling of temperature and density profiles from other biases, we 
compare the mean mass bias from polytropic /3 and extended 
y6-model with spectroscopic temperatures only when no extrap- 
olation of the profiles is needed. Notice that the two models 
are based on projected and three-dimensional spectroscopic-like 
temperature profiles, respectively. Nevertheless, the difference 
between the two types of profiles is small (see Sect. 13. 4t and 
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Table 2. Mean values and standard deviations (errors) of the mass bias derived from the observational-like analysis. These are hsted 
in percent for the three adopted models, the three different fitting ranges, and the three overdensities at which they are evaluated. 
Results are given for the whole sample and for the four 113 classes. 



bias 


fitting range 


rA 


whole sample 


91 


92 


13 


94 






''2500 


-12± 11 


-9 ±6 


-10±9 


-12 ± 10 


-15 ± 14 




[0.1 -0.5]r2oo 


''500 


-27 ± 1 1 


-21 ±7 


-24 ± 10 


-30± 12 


-32± 11 






'200 


-32 ± 12 


-26 ±7 


-30± 13 


-33 ± 11 


-41 ± 12 






''2500 


-14 ±9 


-11 ±6 


-13±8 


-15± 10 


-17 ± 11 




[0.1 -0.75]r2oo 


'500 


-29 ± 1 1 


-23 ±7 


-26 ± 10 


-33 ± 12 


-36 ± 10 




''200 


-34 ± 12 


-25 ±7 


-33 ± 13 


-35 ± 10 


-43 ± 1 1 






''2500 


-16 + 8 


-13±5 


-15±8 


-18±8 


-20 ± 10 




[0.1 - IJrjoo 


'500 


-31 ± 12 


-24 ±7 


-27 ± 10 


-35 ± 13 


-37 ± 10 






''200 


-34± 13 


-25 ±7 


-33 ± 15 


-35 ± 12 


-42 ± 10 






''2500 


-7± 10 


-3 ±8 


-7 ± 8 


-9± 10 


-11 ± 11 




[0.1 -0.5]r2oo 


''500 


-20 ± 23 


-5 ±26 


-16 ±23 


-30 ± 15 


-30± 18 






''200 


-28 ± 42 


-5 ±52 


-30 ± 35 


-35 ± 38 


-42 ±31 






''2500 


-7± 10 


-4± 8 


-6 ±7 


-8± 11 


-10± 11 




[0.1-0.75]r2oo 


''500 


-17+ 19 


-7 ± 15 


-14 ± 15 


-21 ±20 


-28 ± 16 






''200 


-24 ± 26 


-7 ±21 


-22 ± 25 


-29 ± 26 


-41 ± 18 






''2500 


-7± 10 


-4± 8 


-7 ±8 


-9± 11 


-11 ± 11 




[0.1 - l]r2oo 


'500 


-17± 15 


-10± 11 


-13 ± 13 


-21 ± 16 


-24 ± 14 






''200 


-21 ± 20 


-10± 14 


-19± 19 


-23 ±21 


-33± 19 






''2500 


-2± 10 


-2 ±8 


-3 ±8 


-4± 11 


0± 10 




[0.1-0.5]r2oo 


''500 


-13 ±22 


-1 ±27 


-9 ±21 


-21 ± 13 


-20 ± 17 






''200 


-20 ± 42 


1 ±55 


-25 ± 30 


-25 ±41 


-30 ±31 






''2500 


-1 ± 10 


-2 ±7 


-3 ±8 


-2± 12 


1 ± 11 




[0.1-0.75]r2oo 


'500 


-10± 16 


-2± 13 


-7± 13 


-12± 18 


-18± 16 




'200 


-14 ±24 


1 ±20 


-12 ±25 


-17 ±24 


-30± 16 






''2500 


-2± 10 


-2 ±7 


-3 ±9 


-2± 12 


0± 11 




[0.1 - l]r2oo 


''500 


-9± 12 


-5 ±9 


-6± 10 


-12± 14 


-13 ± 12 






'200 


-12± 17 


-3 ± 12 


-11 ± 15 


-11 ± 18 


-23 ± 17 



does not introduce any significant additional bias, in particular 
for the relaxed clusters. 

We therefore focus on the filled symbols in Figs. [T] and |2] 
and the corresponding values in Table |2] It is important to re- 
mark that, as ubiquitously found in simulations, the gas density 
slope in our simulated clusters considerably steepens with ra- 
dius. In addition, also radial temperature profiles show a non- 
trivial radial dependence. Thus, the performance of a given 
model adopted for the surface brightness and temperature fitting 
strongly depends on the ability to model these slope changes. 

From a visual inspection of the fitted profiles we find that, as 
expected, the extended y6-model is extremely accurate in mod- 
eling surface brightness and temperature profiles for any given 
radial fitting range, especially for the most relaxed clusters. The 
accuracy of this model is reflected in the fact that, when no ex- 
trapolation is involved, the bias mean values and standard de- 
viations do not depend on the radial fitting range (compare the 
filled symbols in the different panels in Fig. |2] and the values in 
in Table |2]i. 

The difference between the mean values derived from ex- 
tended /3- and polytropic j6-models (compare filled symbols in 
the bottom panels of Figs. [T] and |2] or the corresponding values 
listed in Table |2]i provides a direct measure of the bias due to 
the inaccurate modeling provided by polytropic j6-model. The 
comparison shows that the polytropic y6-model introduces an ad- 
ditional mass underestimate which is on average ~ 5, 10, and 
15% at r250o, '"500, and r2oo, respectively. 

4.2.2. The extrapolation bias 

The bias originating from extrapolating of the mass profiles be- 
yond the radial range probed by observation can be extremely 



large. While this type of bias can be easily kept under control 
by simply avoiding any mass estimation/comparison beyond a 
given outer radius, it may be in practice the cause of many dis- 
agreements between different mass estimates (i.e. different X- 
ray analyses. X-ray and tensing, etc.). 

The bias introduced by extrapolation of course affects all the 
explored models. However, in order to evaluate its net effect, 
we restrict our discussion to the results obtained from the ex- 
tended y6-model and mass-weighted temperature profiles, since 
spectroscopic temperatures introduce an additional bias that is 
discussed below. The polytropic y6-model, in addition to being 
inaccurate (see Sect. |4.2.T] ), also affected by this latter bias. For 
a given overdensity, the average extrapolation bias can be there- 
fore derived by comparing open with filled symbols in Fig. [3] and 
the corresponding values in Table|2] 

While the extrapolation bias is of course due to the fact that 
best fit functions for data in a given radial range do not provide 
a good description of data outside the radial range, it is inter- 
esting to notice that it causes a systematic underestimate of the 
total mass. For the whole sample, the difference between mean 
mass bias obtained adopting the the fitting range [0.1-l]r2oo (i-C- 
with no extrapolation) and those obtained from the fitting range 
[0.1 - 0.5]r2oo indicates that extrapolation causes an additional 
~ 10% average underestimate at r2oo (see Table [2]i. For a very 
large sample the average bias introduced by extrapolation is thus 
not extremely large. Furthermore, if only the most relaxed clus- 
ters (i.e., the sub-sample ^1) are taken onto account this bias is 
fully negligible (see Table |2|i. These considerations are true for 
average bias values. It is however extremely important to notice 
that standard deviations of mass biases computed from extrap- 
olated values are at least twice as big as those obtained from 
non-extrapolated values, independently on the dynamical state 
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Fig. 2. Same as in Fig.[Tl but for the mass reconstruction method 
adopting the extended yS-model and 3D spectroscopic-Uke tem- 
perature profiles. 



of the cluster (see the errorbars in Fig. |3] and errors in Table |2]i. 
This clearly shows that, for individual objects, extrapolation can 
lead to extremely large mass over/underestimates even for the 
most dynamically relaxed objects. 

4.2.3. The dynamical state bias 

X-ray observations of unrelaxed clusters are expected to yield 
bias mass estimates because both assumptions of spherical sym- 
metry and hydrostatic equilibrium are not valid. 

Here we quantify the mass bias due to the unrelaxed dynam- 
ical state of the cluster by considering: estimates derived from 
the extended yS-model (because of its accuracy, see 14.2.11 ). 3D 
mass-weighted temperature profiles (in order to avoid spectro- 
scopic temperature biases), and mass estimates derived with- 
out extrapolation (to avoid contamination from extrapolation bi- 
ases). Our results are shown by the filled symbols in Fig.[3]and 
since for the case under consideration the results do not depend 



Fig. 3. Same as in Fig. |2l but for results derived using mass- 
weighted temperature profiles. 



on the adopted fitting range we focus our discussion on the val- 
ues shown in the lower panel (i.e. the fitting range [0.1 - l]r2oo, 
see Table|2]for the corresponding values). 

While at r25oo we do not find any appreciable variation of 
the mean mass bias for the four Ylj, classes (see Sect. 13.21 for 
the definition of these classes), at rsoo and raoo we find a clear 
trend with cluster dynamical state: the mass is on average more 
underestimated for disturbed clusters than for the relaxed ones. 
We find, for example, »i(i'exij3T™ ('"200)) = -3% and -23% for 
the qi class (the most relaxed clusters) and q4 class (the most 
disturbed clusters), respectively. 

Furthermore, we find that at all overdensities the individual 
mass bias values are more scattered around the mean for dis- 
turbed objects than for the relaxed ones, as shown by the stan- 
dard deviations in Fig. [3] and errors in Table [2] This is illus- 
trated more specifically in Fig. |4] for the case of fitting range 
[0.1-l]r2ooandA = 500. 

These results are particularly relevant because they clearly 
show that mass underestimates in clusters that are identified as 
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Fig. 4. Mass bias as a function of Ilj at rsoo with estimated 
masses derived adopting the extended j6-model and 3D mass- 
weighted temperature profiles. The solid line shows the best fit 
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>(r5oo) = -0.256 - 0.022 x n3(r5oo). 



relaxed according to the power ratio method are on average small 
(~ 2 - 5%), as opposed to the very large bias that affects total 
mass determinations in disturbed systems . Finally, we remar k 
that, in agreement with other studies (e.g. iJeltema et al.ll2008l) . 
we do not find any correlation between 1X3 and true total mass. 

4.2.4. The spectroscopic temperature bias 

The difference between emission-weighted, mass-weighted, and 
spectroscopic ICM temperatures has been address ed by a num- 
ber of authors (e.g. [Mathiesen & Evrard 20 0][ iGardini et al] 
2004^, 'Mazzotta et al.| |2004£ iRasia et alj|2005l; I Vikhlininl 120061: 



Vald arnini 2006). Special emphasis has been given to global 
temperatures because of their impact on the scaling relations 
such as the luminosity-temperature and mass-temperature rela- 
tions. Recent studies indicate that rn,w(< ^a) < rsi(< ?"a) < 
Jew(< r^ ). where usually A = 500 (e.g. iNagai et al.1 l2007al: 
ICiotti & Pellegrini 2008). This relation is also valid for our sim- 
ulated sample. By averaging over the whole sample we find the 
ratios: r„,w(< '•500) : T,i{< rjoo) : rew(< r^oo) = 1 : 1-13 : 1.27. 
In addition we computed average temperatures from the same 
sample but without applying the masking procedure described 
in Sect, [ll] and find: Tm^{< rsoo) : Ts\{< r^oo) : rew(< rsoo) = 
1 : 1.08 : 1.21. Considering that mass-weighted global temper- 
atures are fully unaffected by masking, these results show that 
both emission-weighted and spectroscopic-like global tempera- 
tures are on average biased high by ~ 5% if cold ICM clumps 
are not masked out appropriately. 

In the following we focus on the difference between Tnm 
and Tsi since emission weight ed tempe r atures are not relevant 
in our analysis. As noted by Na gai et all (l2007al) global spectro- 
scopic temperatures within large apertures are higher that mass- 
weighted ones because they are dominated by the inner, hotter 
cluster region. It is important to notice that accurate total mass 
determination methods such as the one adopted in this paper 
rely on spatially resolved temperature determinations and not on 
global temperatures. Therefore it is crucial to compare T^^ and 
Tsi radial profiles and how any difference between the two prop- 
agates into total mass estimates. 

We compared spectroscopic like and mass-weighted temper- 
ature profiles and find, in agreement with Ameglio et al. ( 200^, 
that Tjawir) > T^iir). This inequality, and the fact that the op- 
posite is true when considering global temperatures, can be un- 



derstood as follows. Global temperatures are computed by av- 
eraging within large cluster-centric distances whereas profiles 
are computed by averaging within radial shells. The density and 
temperature distributions in the two cases are dramatically dif- 
ferent. In the former case the distribution is dominated by the 
central, hot and dense gas, while in a given radial shell the ICM 
is obviously distributed in a much narrower range of densities 
and temperatures. The ICM in our simulations is locally mul- 
tiphase and made of a dominant component and mostly lower 
temperature particles which do not substantially contribute to 
mass-weighted temperature, but bias spectroscopic temperatures 
low. A simple illustration of the spec troscopic tempe rature bias 
for a two-phase ICM can be found in lJia et alj ( l2008h . 

It is very relevant to remark that in the computation of both 
Trnwit") and Ts\{r) all cold clumps are excluded. This implies that 
local spectroscopic-like temperatures can underestimate mass- 
weighted ones also because of cold particles which are not part 
of clumps. Interestingly, we find that our mass reconstruction 
method depends very weakly on cold clumps removal. Total 
masses are in fact derived by fitting brightness and temperature 
profiles with smooth functions of radius, which are in almost all 
the cases insensitive to the irregularities (spikes or bumps in the 
brightness profiles and dips or depressions in the temperature 
profiles) arising from cold clumps when no masking is applied. 

In relaxed clusters T^^/ and T^i agree very well at small 
and intermediate radii. The agreement becomes gradually worse 
with increasingly disturbed dynamical state (as measured by 
the power-ratio method). This reflects the more complex ther- 
mal structure of disturbed clusters, where the fraction of cold 
gas particles stripped from cold sub-structures is relatively high. 
Furthermore, we find that the relative difference between Tj„„{r) 
and Tsi{r) increases with radius. This is expected because of the 
weighting scheme used in the computation of Tsi and the pres- 
ence of cold infal ling gas in the clus ter outskirts. A similar find- 
ing is reported inl Rasia et al.l (l2006h . Finally, the relative differ- 
ence between and T^i is higher in more massive clusters than 
in low mass systems be cause in massive objects the spread i n 
temperature is larger (e.g. lValdarninil2006l:lAmeglio et alj2007h . 

The underestimate of the true ICM temperature profile in- 
troduced by spectroscopic measurement leads to total mass un- 
derestimates because hydrostatic masses (Eq. flAi ) depend on 
both temperature and its logarithmic derivative. In order to dis- 
entangle this effect from the other observational biases discussed 
in the previous subsections, we consider results obtained from 
the two extended yS-models and compare results obtained from 
the one adopting spectroscopic-like profiles with those derived 
using mass-weighted temperature profiles. The net mass bias in- 
troduced by spectroscopic temperatures can be readily estimated 
by computing the difference between the mean bias values listed 
in Table |2] However, a clearer determination can be performed 
by defining the spectroscopic bias: 



M"" ,„(< r) - M'" . 



(< r) 



.{<r) 



(16) 



where M''"' ,„(< r) and M*"" ,„(< r) are total masses esti- 

mated through the extended yS-model with spectroscopic-like and 
mass-weighted temperatures profiles, respectively. Both are de- 
rived adopting the the fitting range [0.1 - l]r2oo in order to avoid 
any contamination introduced by the extrapolation. 

In Table |3] we list the spectroscopic bias mean values and 
standard deviations computed for the whole sample and for the 
four 113 classes separately. These values are also shown in Fig. 
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fitting range: [0.1-1]r2oo 



q2 

113 (r^ ) classes 

Fig. 5. Mean (points) and standard deviation (errorbars) of the 
spectroscopic bias i'.spec('") (Eq. ( fT6l )) at rom (stars), rsoo (di- 
amonds), and r25oo (triangles). The results are shown for the 
whole sample (all) and when the four TI3 classes are used as 
sub-samples. The results for different overdensities are shifted 
horizontally to improve clarity. 
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Fig. 6. Mass bias from the extended yS-modeling versus total 
true mass at A = 500 for the radial fitting range [0.1 - l]r2oo- 
The estimated total mass is derived using mass-weighted tem- 
perature profiles. The solid line shows the best fit function 
^..«/?,7-s('"50o) = 0.78 - 0.061og[M(< rsoo)//^' Mq]. 



|5] for the sake of clarity. The spectroscopic bias is systematic 
and of the order of -10% if the whole sample is considered. 
However, if the same comparison is performed for the four Tl^ 
classes separately one finds that it is larger (in absolute value) 
in disturbed clusters than in relaxed ones. The mean spectro- 
scopic bias is also more important at larger cluster-centric dis- 
tances than at small ones. These results are of course directly 
related to the relative difference between the T^v, and T^i pro- 
files discussed above. 

To reiterate, a crucial point of out findings is that also after 
the removal of all the resolved ICM cold clumps in the simula- 
tions, spectroscopic temperature profiles underestimate the true 
ones, and so the total mass, because of a diffuse distribution of 
cold gas particles. 

4.2.5. Mass dependence 

Our results concerning the influence of the various observational 
biases allows us to focus on some relevant cases when exploring 
the dependence of the mass bias on the total, true mass. We can 
in fact avoid any contamination from biases such as the model- 
ing bias or the extrapolation bias by considering results provided 
by extended j6-models and without extrapolation of the mass pro- 
files. 

As expected, no mass dependence is found for mass biases 
derived when using the extended j6-model, as shown in Fig.|6]in 
the case of fitting range [0.1 - l]r2oo and A = 500. The indepen- 
dence of bias on total mass is found for all fitting ranges and all 
overdensities. This is, however, only found when mass-weighted 
temperature profiles are used. Spectroscopic-like temperatures 
introduce a mass dependence: the average bias is more nega- 
tive for massive systems, as shown in Fig. [T] for a specific case. 
This mass dependence arises from the weighting scheme used 
in the computation of spectroscopic-like temperatures, which is 
sensitive the cooler part of the gas distribution, and the scale 
dependency introduced by cooling on the amount of this cold 
component versus the cluster mass (see Sect. 14.2.41 ). 

In order to investigate the dependence of the spectroscopic 
bias on true mass and at the same time on dynamical state, we 
compute mean values and standard deviations by considering 
objects in a given (true) mass interval and 1X3 class. In order to 



obtain statistically significant values (i.e. a sufficient number of 
object in each class), we must adopt a coarse mass binning. After 
experimenting with different types of mass binning we find that 
the discussion can be greatly simplified by reporting results for 
only two sub-samples, which are obtained by ordering the sam- 
ple by true cluster mass and halving it. In the following these 
two sub-samples are referred to as lowM- and /j/g/jM-sample. 

For each of these sub-samples we construct the four Yl^ 
classes separately and compute mean and standard deviation of 
the spectroscopic bias as done above for the whole sample. The 
results are reported in Table [3l together with those obtained for 
the whole sample, and shown in Fig. [8] At r25oo we do not find 
any significant dependence on mass and therefore the average 
biases reported in Sect. 14.2.41 are valid at any mass. On the other 
hand, we find a strong mass dependence at rsoo and r2oo and 
that its strength increases with cluster-centric distance. For the 
/owM-sample the average spectroscopic bias does not signifi- 
cantly vary with overdensity, as opposed to the behavior found 
for the /i/^/iM- sample. Because of the findings reported in Sect. 
14.2.31 relaxed clusters {q\ and ^2 sub-samples) are of particular 
relevance. The ones in the ZowM-sample are very mildly affected 
by the spectroscopic bias at all cluster-centric distances, whereas 
for those in the /i/g/jM-sample the mean spectroscopic bias is 
10% at rsoo and ~ -15% at r2oo- 

4.2.6. The ideal case 

After having examined the various biases affecting the total mass 
determination, we discuss here the ideal conditions under which 
the most unbiased mass determination is possible. From our re- 
sults it is clear that in order to obtain unbiased mass estimates 
one must select very relaxed clusters (the q\ sub-sample), adopt 
a model that can accurately fit the ICM density and temperature 
profiles (the extended j6-model, here), and avoid extrapolation 
(notice, however, that the extrapolation bias is on average neg- 
ligible for relaxed clusters. See Sect. |4.2.2t . In addition let us 
assume that the true thermal structure of the ICM (i.e. rn,„(r)) 
can be measured. In this case the bias does not depend on the true 
mass of the system and, at all radii, masses are on average very 
well reconstructed (see Table[2]and Fig.O. Notice that there is 
a residual level of underestimation which might be originating 
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Table 3. Mean values and standard deviations (errors) of the spectroscopic bias /7spec('"A) (Eq. (fT6Tl). Values are in percent and 
listed for three different true mass ranges used in their derivation: all masses (the whole sample), /owM-sample (clusters with 
M(< r2oo) < 1-8 X IO'^^/i^'Mq), and highM-sample (clusters with M(< r2oo) > 1.8 x lO'^^/r'Mo). In each case results are given, at 
three radii, for the whole sample and for the four TI3 classes. 



mass range 


rA 


whole sample 


91 


12 


93 




94 




''2500 


-6± 


6 


-2 + 2 


-4 ±4 


-6± 


5 


-11 ±7 


all masses 


''500 


-9± 


7 


-5 + 4 


-8 ±5 


-12 d 


:7 


-13 ±8 




''200 


-11 ± 


10 


-7 + 6 


-9 ±8 


-14 + 


12 


-14 ± 10 




''2500 


-4± 


4 


-1 ± 1 


-3 ±3 


-6 + 


4 


-9 ±5 


ZowM-sample 
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-7 ± 
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-4 ±2 


-5 ±3 


-9± 
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-9 ±4 
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-3 ±2 


-4 ±4 


-7± 


5 


-9 ±5 
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-7 ± 
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-2± 1 


-4 ±4 


-7± 
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-12±9 


highM-sample 
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-12 ±8 
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-11 ±6 


-14 d 


:8 


-16±9 
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-12±6 


-17 ±7 


-19 ± 


14 


-18 ± 11 
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Fig. 7. Same as Fig.|6] but with the estimated total mass derived 
using spectroscopic-like temperature profiles. The solid line 
shows the best fit function b^^^ 7--'i)('"5oo) = 2.79 - 0.21 log[M(< 
r5oo)//r'M0]. 
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Fig. 8. Same as Fig. |5] but for the lowM- (filled symbols) and 
/i/g/jM-samples (open symbols), separately. 



from the still imperfect gas density and temperature modeling or 
by the departure of the gas from the hydrostatic equilibrium. 

It is of paramount importance, however, to take into account 
that the true ICM temperature is not accessible through X-ray 
observations and that we must deal with spectroscopic temper- 
ature measurements. As shown in Sect. 14.2.41 the most relaxed 
clusters are also the less affected by the spectroscopic bias (see 
Table[3]and Fig.lSj. As we have shown in Sect. 14.2.31 the spec- 



troscopic bias depends on cluster mass. It is negligible in low 
mass clusters, but important for massive ones especially at large 
radii. 



5. Total mass determination from the 
three-dimensional gas distribution 

As shown in Sect. |4] mass determinations from X-ray observa- 
tions suffer from various biases. Of course, in addition to these 
also the assumption of hydrostatic equilibrium (i.e. Eq. ( fT4b ) 
plays a crucial role in the total mass reconstruction. Moreover, 
the derivation of Eq. ( fT4l ) implicitly assumes also spherical sym- 
metry. Whereas the validity of this latter assumption can be fairly 
well tested though an X-ray morphological analysis, the robust- 
ness of the assumption of hydrostatic equilibrium can not be di- 
rectly investigated observationally. 

Here, in order to investigate the validity of both assumptions 
we avoid any bias due to the X-ray reconstruction method and 
make use of the full three-dimensional gas distribution in the 
simulation. For consistency with the analysis presented in Sec.|4] 
we consider only objects with rsi(< r2oo) ^ 2 keV. 

5.1. Mass estimators 

We introduce the hydrostatic equilibrium mass estimate, Mhe 
by assuming spherical symmetry and rewriting Eq. ( fT4b as 

Mhe{< r) = -— < V,P,(r) >, (17) 

where Pg is the gas pressure. As a means to avoid any system- 
atic effect generated by estimating gradients from averaged pro- 
files and to achieve maximum accuracy, we evaluate the term 

< V,Pi;(r) > directly from the gas particle distribution. To this 
end, we introduce a spherical shell at the test radius r and a cor- 
responding set of 40 X 40 grid points with angular coordinates 
uniformly spaced in cos9, (p. At the grid point Xg,- the pressure 
gradient is expressed as 

VPgiXgr) = (r - 1) X "^'■"'■'^^(l^' - ^grl hg), (18) 

where m, is the specific particle internal energy and j - 5/3 
. The radial component is then extracted by transforming 
the gradient vector according to the grid coordinates and finally 

< V, Pg(r) > is obtained by averaging ( fTSl ) over the set of grid 
points. We checked the robustness of the estimated average gra- 
dients < Vff g(r) > with a two-fold increase in the number of 
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grid points (i.e. 80 x 80 ) and found the results to be unaffected 
by this choice. 

This procedure allow us to define a mass estimate Mhe 
which is influenced by the cluster dynamical state, but in con- 
trast to the mass determinations discussed in Sect. |4] it is not 
affected by any bias originating from the gas density and tem- 
perature profiles reconstruction and modeling. 

Similarly to Eq. ( fTSl l we define the mass bias: 



Mhe{< r) - M(< r) 
M{< r) 



(19) 



In addition, w e also consider the hydrodynamic al equilibrium 
mass estimate dRasia et al ] l2004HKav et"aLll2004 : 



Mde(< r) 



Mhe(< r) + MA< r) 
I d In pg 

g'' 



Mhe{< r)-^ 



d Inal 
+ ^^+2/?v(r) ,(20) 



din r d In r 



where cr,- and cr, are the radial and tangential gas velocity dis- 
persion, respectively, and y6,,(r) = 1 - crj llcr^ is the gas velocity 
anisotropy parameter The estimate of the mass term M,j{< r) 
is subject to uncertainties because of the irregularities which are 
present in the gas velocity dispersion profiles of the simulated 
clusters. In order to properly estimate the radial derivatives in 
Eq. |20] we have therefore regularized the me asured profiles by 
applying to them a Savitzky-Golay filter (e.g. iPress et alj[r992l 
Sect. 14.2). This smoothing procedure is effective in removing 
the small-scale noise and yields mass estimates which are much 
more regular against the radial dependency. As for Mhe we de- 
fine the mass bias: 



boEir) = 



Mde{< r) - M(< r) 
M{< r) 



(21) 



Once the mass profiles Mhe{< f) and Moe{< f) for the ~ 100 
objects in our sample are computed, we evaluate the mass biases 
bHsir) and bosir) at ^2500, '"500, and r^oo- Finally, mean values 
and standard deviations of the resulting mass bias distributions 
are computed. 

Notic e that Eq. (pO| is derived f rom the Boltzmann equa- 
tion (e.g. iBinnev & Tremain^ Il987l Chapter 4) by assuming 
steady-state hydrodynamic equilibrium and taking into account 
the terms from the gas isotropic pressure and anisotropic veloc- 
ity dispersion. If the latter are neglected (i.e. it is assumed that 
the velocity dispersion of the gas is much smaller that its tem- 
perature), then Eq. ( l20l i yields Eq. ( fTTI ). In the derivation it is 
assumed, of course, that the system is spherically symmetric. 

The steady-state assumption implies no net radial streaming 
motions (no mean laminar flow, i.e. < v, >= 0, which implies 
crj =< vj. >) and that all time derivatives can be set to zero. 
Therefore, the applicability of Eq. ( l20l i depends on the degree of 
spherical symmetry of the system and on the condition < v, >= 
0, while that of Eq. ( fTTT i additionally depends on the condition 
Ma-{< r)IMt^{< r) «: 1. Hence, the hydrodynamical equilibrium 
equation naturally provides a more accurate mass estimator than 
the hydrostatic equilibrium. 

The accuracy of spherical symmetry can be assessed by mea- 
suring the mean tangential acceleration at the surface of the 
spherical shell of radius r. To this end, in correspondence of each 
of the grid points at which we evaluate the pressure gradient in 
Eq. ( fTSl l. we also evaluate the gravitation acceleration Ug of test 
particle. We then define a mean torque parameter Tg by averag- 
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Fig. 9. Mean torque parameter (which measures the spheric- 
ity of the gravitational potential) as a function of < TIt, >pianes 
(which measures the sphericity of the gas distribution) at r25oo. 
For each cluster the value < 113 >pianes is the mean of the three 
113 values computed from three orthogonal surface brightness 
maps. 



ing the acceleration components Ug + and over the set of 
grid points: 



ii:, 



7 1 



(22) 



The constraint Tq <«; 1 can be used as a condition to validate the 
spherical symmetry approximation. Accordingly, the quadrupole 
moment terms can be neglected in the radial component of the 
total potential gradient and to the lowest order we have: 



V,<l)(r) 



GM{< r) 



Tq « 1. 



(23) 



It is interesting to notice that the degree of spherical symme- 
try measured directly from the full 3D information provided by 
the simulations is highly correlated with that measured from X- 
ray brightness maps. As shown in Fig.|9]for the specific case of 
A = 2500 the torque parameter Tq and < 113 >pianes are corre- 
lated. Here < 113 >pianes is the mean of the three 113 values com- 
puted from three orthogonal surface brightness maps (see Sect. 
13. 2t . From a Spearman rank correlation coefficient analysis we 
find that the correlation is significant at more than the 95% level 
at r25oo, '■500, and r2oo- Clusters with small values of torque pa- 
rameter Tq, i.e. whose potential has a high degree of sphericity, 
are therefore identified as spherical systems when selected by 
means of the power ratio method. 

The level of radial streaming motions at the radius r is 
directly quantified through the radial Mach number < v,. > 
/cj, where the gas sound speed Cj, is computed from Cj(r) = 
5fcrg(r)/(3//OTp) 

It follows that the applicability of Eqs. ( fTTj i and ( [2Q] i to mea- 
sure the cluster gravitating mass within the test radius r is then 
limited to those systems for which at the radius r the conditions 
T„ «; 1 and | < > /c,! «; 1 are satisfied. 



5.2. Results 

The dependence of the mass biases bnE and boE on the mean 
torque parameter Tg and the radial Mach number < v, > /c, are 
shown in Figs. [TOl and fTTI respectively. The figures show that. 
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for a given overdensity, we measure large mass biases (both pos- 
itive and negative) for clusters with large Tg and | < v,. > /cj|. 
As we move from large to small values of the parameters (i.e. 
to an increasing degree of sphericity of the cluster or to small 
radial streaming velocities) the scatter of the measured biases 
substantially decreases. These results demonstrate that, in addi- 
tion to sphericity, the radial Mach number < v,. > /c^ is a key 
parameter in the study of the reliability of the hydrostatic and 
hydrodynamical mass estimators. 

It is important to notice that the parameters Tg and < v,. > /c., 
are not uncorrelated. We find in fact that objects with large 
also tend to have large | < v, > /c,!, as shown in Fig.[T2]for A = 
2500. The correlation between these parameters is found to be 
significant at more than the 95% level at all the three considered 
overdensities. 

In addition to the considerable decrease of the scatter with 
decreasing values of Tg and | < v,. > /c^l, a visual inspection 
of Figs.fTOlandfTTIindicates that on average the mass estimators 
Mhe is biased low and that masses computed through the more 
accurate estimator Moe are on average well recovered. This is 
quantitatively shown by mean values derived from the whole 
sample which are listed in Table |4] 

It is crucial to remark that average biases derived from the 
hydrostatic mass estimator Mhe are in agreement with those ob- 
tained from the X-ray mass reconstruction procedure presented 
in Sect. |4] In fact, a comparison between the values in Table 
ID and those listed in Table |2] for extended yS-model with mass- 
weighted temperature profiles (since these follow the virial re- 
lationship more accurately that spectroscopic ones) shows that 
in the two cases the true mass is similarly underestimated. 
Furthermore, we find that, as for the mass bias derived using 
mass-weighted temperature in Sect. |4] both bnE and boE do not 
depend on cluster true mass. 

Most important, the mean values of the bnE and bos distri- 
butions show that if in addition to the thermal pressure also the 
ICM non-thermal pressure component (i.e. the gas anisotropic 
velocity dispersion) is taken into account, the total mass of the 
system is on average better reconstructed (see also iRasia et al.l 
l2004tlKavetani2004 . 



Table 4. Mean values and standard deviations (errors) of the 
mass biases and boE derived from the whole sample and for 
the sub-sample with Tg < 0.15 and | < iv > /c.s| < 0.1 . Values 
are given in percent at three dififerent overdensities. 
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In Oder to investigate the accuracy of the mass estimators 
Mhe and Mde more precisely we proceed as follows. At a given 
overdensity, we create sub-samples of clusters with Tg and | < 
Vr > /cs\ smaller that some given threshold values and in each 
case we compute mean values and standard deviations of bHE 
and boE ■ It is important to find a compromise between small 
threshold values for Tg and | < v, > /c j (in order to guarantee 
as much as possible the applicability of the mass estimators) and 
a sub-sample with enough objects so as to allow a meaningful 



Fig. 10. Mass biases bnE (open stars) and boE (filled triangles) 
as a function of mean torque parameter Tg at r25oo, I'soo, and r2oo- 
In each panel we show the mean values derived from the whole 
sample for bnE (dashed line) and boE (dot-dashed line), and the 
mean values (plotted only up to 0. 1) derived from the sub-sample 
with Tg < 0.15 and \ < Vr > /cs\ < 0.1 for bnE (short-dashed 
line) and boE (solid line). In each panel are also reported these 
mean values along with the standard deviations. Notice that in 
order to focus on the bulk of the data we have omitted the results 
for a few clusters with very large bias and large Tg. 



distribution fitting. We experimented with various cuts and find 
that the optimal choice is to select clusters with Tg < 0.15 and 
I < V,. > /cjI < 0.1, which yields sub-samples of at least ~ 40 
objects at each overdensity. Results derived from this sub-sample 
are listed in Table |4] We remark that the results are not sensitive 
to exact choice of the threshold values. As expected the standard 
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Fig. 11. Mass biases as a function of radial Mach number < > 
/c,. Symbols, lines, and values are as explained in Fig. [10] Also 
here few clusters with very large bias and large | < v, > /cj| are 
omitted from the plot to improve clarity. 



deviations are much smaller than those derived for the whole 
sample. Standard deviation for the sample with > 0.15 and 
I < V,- > /CjI > 0.1 are ~ 2 - 3 times larger that those obtained 
for the sample with < 0.15 and | < > /c^l < 0.1. 

The analysis of simulated objects for which the applicability 
criteria of Eqs. ( fTTI i and ( l20l i are most valid allows us, therefore, 
to draw more robust conclusions on the reliability of the mass 
estimators Mhe and Moe- Hydrostatic masses are biased low 
by ~ 5% at r25oo and by ~ 10% at rsoo and r2oo, whereas the 
hydrodynamical equilibrium mass estimator is unbiased at any 
overdensity. 
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Fig. 12. Mean torque parameter as a function of Mach number 
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6. Discussion and conclusions 

This work is aimed at elucidating the reliability of X-ray total 
mass estimates in clusters of galaxies using N-body/SPH simula- 
tion of a large sample of clusters. The physical modeling the gas 
includes radiative cooling, star formation, supernovae heating, 
and metal enrichment. The large number of simulated clusters 
enables us to derive very robust conclusions through a statisti- 
cal analysis of the sample. The total mass is recovered adopt- 
ing an observational-like approach and compared with the true 
mass in the simulations. Surface brightness and temperature pro- 
files, that we generate from the simulations, are used to estimate 
the cluster mass at different overdensities (^2500, rgoo, and r2oo) 
by means of the hydrostatic equilibrium equation. We explore 
various models and conditions under which the mass is recon- 
structed in order to entangle different mass biases. In addition, 
a power ratio analysis of the surface brightness maps allows us 
to assess the dependence of the mass bias on cluster dynamical 
state. Moreover, we perform a study on the reliability of hydro- 
static and hydrodynamical equilibrium mass estimates using the 
full three-dimensional gas distribution in the simulation. 
In the following we list and discuss our main findings. 

1. Our analysis shows that it is very important to use ana- 
lytical models with a large amount of parametric freedom 
when modeling the shape of the ICM temperature and sur- 
face brightness radial profiles. A model with a low degree of 
sophistication such as the polytropic y6-model can introduce 
a very large modeling bias. Compared to the more sophis- 
ticated extended /?-model, which is found to be extremely 
accurate in following the slope changes of the gas profiles, it 
additionally leads to average mass underestimates of the or- 
der of ~ 5, 10, and 15% at r25oo, '"500, and r2oo, respectively. 

2. The bias originating from extrapolating of the mass profiles 
beyond the radial range probed by observation can be ex- 
tremely large. We find that the underestimate from extrapo- 
lation alone is of the order of ~ 10% at r2oo, and lower at 
smaller cluster-centric distances, when considering an aver- 
age over the whole sample. However, for individual objects, 
the extrapolation bias can be as large as ~ 50%. 

3. The unrelaxed dynamical state of a cluster can also lead 
to mass underestimates. The total mass is on average bi- 
ased lower for disturbed clusters than for the relaxed ones. 
Furthermore, we find that the bias values are much more 
scattered around the mean for disturbed objects than for 
the relaxed ones. If mass-weighted temperature profiles are 
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adopted, the mean mass bias is at most ~ -5% at all cluster- 
centric distances for relaxed clusters, but can be as large as 
~ -20% for the most disturbed ones. 

Our results a r e in ex cellent agreement with the findings of 
iJeltema et alj (|2008|) . as shown by the comparison of Fig. 
|4] in this work and Fig. 9 (top panel) in their paper, where 
the correlation between the same quantities are shown. This 
agreement i s of particular impor tance, considering that the 
analysis by iJeltema et al.l ( l2008l) is based on simulations 
performed with the adaptive mesh refinement code Enzo. 
Our findings are also in good agreement with the results of 
lKavetal.l(l2007l) . 

Mass estimates derived using spectroscopic temperatures are 
lower that those derived from mass-weighted temperatures 
(i.e. the true gas temperatures). We find that the spectro- 
scopic temperature bias, i.e. the bias originating from spec- 
troscopic temperatures alone, is of the order of -10% for 
the whole numerical sample. A similar value is found by 
lKavetal.l(l2007l) . 

Mass underestimates derived adopting spectroscopic-like 
temperature profiles (7 and 17% on average at A = 2500 
and 500 for the wh ole sample) are in g ood agreement with 
the values found bv lNagai et alj(l2007a ) (12 and 16% are the 
corresponding values), whose results are based on mock X- 
ray observations derived from simulations performed with 
an Eulerian code. We also find very good agreement by 
performing the same comparison for relaxed and unrelaxed 
clusters separately. We n otice that total mass biases derived 
from the simulations by iNagai et al] (l2007al) and adopting 
mass-weighted profiles are on a verage -7% at rjoo for re- 
laxed clusters ( iLau et al.l |2007|) . The corresponding value 
that we find from our Lagrangian simulations is -5%. 
Considering the different nature of the numerical codes and 
the different implementation of the various physical pro- 
cesses, the agreement is extremely relevant. 
Our re sults are in tension with the findings of iRasia et alj 
(120061) . who find, on average, much stronger biases. We no- 
tice, ho wever, that a fair co mparison is not possible since the 
work bv lRasia et alj ( l2006l) focussed only on 5 clusters. 
Our analysis also shows that the spectroscopic bias depends 
on dynamical state. We find that the spectroscopic bias is 
rather small (-2, -5, and -7% at A = 2500, 500 and 200) 
for the most relaxed clusters and of the order of -12% for 
the most disturbed ones. 

While the mass bias derived from mass-weighted temper- 
ature profiles does not depend on true cluster mass, the 
one derived from spectroscopic-like temperature exhibits a 
strong mass dependence. For clusters with M(< ^200) < 
1.8 X IO'^^/i^'Mq we find that, at all radii, the spectroscopic 
bias contributes on average to the total bias with less than 
7%, and less than 4% if one considers only the most relaxed 
clusters. On the other hand, for clusters with M{< rooo) > 
1.8x lO''* /i^'Mq we find that the spectroscopic bias is larger: 
e.g. -7% and -12% for the relaxed clusters at rsoo and r2oo, 
respectively. 

Even in the ideal case, i.e. when the mass of relaxed clusters 
is estimated without extrapolation and adopting a very ac- 
curate model for the ICM density and temperature profiles, 
total masses are affected by the spectroscopic temperatures. 
In this case the mean total mass bias is ~ (-4, -3), (-12,-9), 
and (-15,-3)% at r25oo,'"500, and r2oo for the (most, least) 
massive clusters in our simulated sample. 
Even if we assume that the true (mass-weighted) ICM tem- 
perature can be reconstructed, masses are biased low (by -2, 



-9, and -12% at A = 2500, 500, and 200, respectively, when 
taking the average over the whole sample). 

8. The latter finding prompted us to investigate the possible 
violation of the hydrostatic equilibrium by using the full, 
three-dimensional information provided by the simulations. 
In agreement with the results from our observational-like 
analysis, we find that the hydrostatic equilibrium assumption 
yields masses underestimated by~ 10-15%. This implies 
that the origin of the bias is not originated by the X-ray re- 
construct ion method . The s ame level of mass bias has been 
found bv lRasia et all (|2004 and lBurns et"aLl (|2008|) . 

9. In order to elucidate the origin of this bias we compute 
masses using both hydrostatic and dynamical equilibrium 
mass estimators. The dynamical equilibrium takes into ac- 
count for both thermal and non-thermal pressure of the ICM. 
We find that the mass estimated through dynamical equilib- 
rium are on average well estimated. T he mean bias pro files 
of the "feedback" clusters simulated bv lKav et al.l (|2004|) are 
in very good agreement with our average values at all over- 
densities. 

10. We explore the conditions of apphcability of both estima- 
tors, i.e. spherical symmetry (torque parameter <K 1) 
and small radial streaming velocities (radial Mach number 
I < \v > /C.5I « 1). We find that the biases bHE and boE 
do not depend on cluster true mass and that their scatter de- 
creases with decreasing Tq and | < v,. > jc^Y 

11. For clusters with small values of Tq < 0.15 and | < \v > 
/cjI < 0.1 we find that the average mass underestimate found 
from the hydrostatic equilibrium estimator (5% at r25oo and 
10% at rsoo and r2oo) is extremely well corrected by adopting 
the dynamical equilibrium estimator 

Our results of course depend on the physics included in the 
simulation. Our modeling of the gas physical processes is in- 
complete and does not include, for example, energy feedback 
from active galactic nuclei (AGN). The inclusion of more real- 
istic physics is particularly relevant in the cluster cores, since it 
is supposed to provide additional heating to the gas which could 
help to solve the cooling flow problem as well as the over cool- 
ingof baryons actually present in current simulations (e.g. lVoill 
120051) . 

The effect on the thermal status of the ICM of energy feed- 
back from the AGN are however not expected to modify in a 
significant way the ICM properties outside the cluster cores. The 
validity of this assumption is justified from the radial behavior of 
the measured temperature profiles, which are in good agreement 
at r > 0.1 X r20() with the present simulations (|yaldarninil l200"6l) 
and the ones of other authors (e.g. lNagai et al.ll2b07bl) . 

Importantly, we find that, in addition to cold gas clumps, also 
the diffuse cold gas component which is left after the removal of 
all resolved clumps substantially biases spectroscopic tempera- 
tures low. The good agreement of our temperature profiles with 
those found by other authors using different codes suggests that 
the amount of the cool gas component present in our simulations 
is not affected by numerical issues. 

There is another issue, which is connected to the gas phys- 
ical modeling in the simulations, that could be potentially rel- 
evant for our analysis. In the runs performed here the artificial 
viscosity is treate d according to the standard SPH formulation 
(lMonaghanll2005l). which is a nume rical scheme comparatively 
VISCOUS orris & Monaghanlll997h and inadequate to follow 
the development of fluid turbulence. It has been shown that the 
level of kinetic energy in random gas motion could be as high as 
~ 30% of the gas thermal energy dPolag et al.ll2005l:IVazza et~al1 
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20061). when the numerical viscosity scheme is generalized as in 
Morris & MonaghanI d 19971) to properly treat fluid turbulence. 

It is worth noticing that very high levels of random mo- 
tions are inconsistent with the upper limits recently derived 
from the comparison of X-ray and weak lensing mass estimates 
jMahdavi et alj|2008l: Izhang et al.ll2008l) . Direct measurements 
of the various sources of non-thermal pressure (gas motions, 
cosmic rays, magnetic fields, etc.) are extremely diflicult, and 
the comparison between X-ray and lensing masses presently 
provides the most significa nt constraints on th e level of non- 
thermal pressure in the ICM. IZhang eTal] (|2008|) obtain a ratio of 
1 .09 ± 0.08 between the weak lensing and X-ray mass estimates 
extracte d at rspp from a sampl e of 19 clusters. At the same over- 
density, jMahdavi et al. I l20^ find that the ratio between X-ray 
and lensing is 0.78 + 0.09 (0.85 + 0.10 after correction for excess 
structure along the line of sight) for a sample of 18 clusters. The 
fairly good agreement between these recent observational results 
and the values found in the present analysis and other similar 
studies therefore implies a low level of turbulence present in the 
ICM, thereby suggesting that the gas physics outside the cluster 
cores is well approximated by the simulations presented here. 

While we have shown that at present both theoretical models 
and observations seem to indicate a deviation from the hydro- 
static equilibrium of the order of ~ 10%, further investigations 
are needed in order to draw robust conclusions on this impor- 
tant issue. Future measurements of X-ray and lensing masses for 
large samples and measurements of the ICM velocity structure 
will provide valuable information on the validity of the hydro- 
static equilibrium and therefore on the reliability of X-ray clus- 
ters as cosmological probes. 

Moreover, these measurements will also be of particular 
relevance for constraining gas physics in galaxy clusters. The 
amount of physical viscosity is in particular a key parameter 
which quenches the level of turbulence present in the ICM. 
Physical viscosity has been incor porated in hyd rodynamic sim- 
ulations of galaxy clusters by iSiiacki & Springel (2006) and it 
was shown that even a modest amount of physical viscosity has 
significant consequences on ICM properties. Therefore X-ray 
mass estimates, extracted from a statistically meaningful sample 
of hydrodynamical SPH simulations which include physical vis- 
cosity, could be profitably used to indirectly measure the level of 
ICM viscosity when contrasted against X-ray and weak lensing 
mass measurements. 

These constraints will also likely have a significant impact 
on those scenarios in which the cooling flow problem is solved 
by providing additional heating to the gas through energy dissi- 
pation. 
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